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ABSTRACT 


When the light from a distant object passes very near to a foreground galaxy or cluster, gravitational 
lensing can cause it to appear as multiple images on the sky [1]. If the source is variable, it can be used to 
constrain the cosmic expansion rate [2] and dark energy models [3]. Achieving these cosmological goals 
requires many lensed transients with precise time delay measurements [4]. Lensed SN are attractive 
for this purpose because they have relatively simple photometric behavior, with well-understood light 
curve shapes and colours—in contrast to the stochastic variation of quasars. Here we report the 
discovery of a multiply-imaged supernova, AT 2016jka (“SN Requiem”). It appeared in an evolved 
galaxy at z = 1.95, gravitationally lensed by a foreground galaxy cluster [5]. It is likely a Type Ia 
supernova—the explosion of a low-mass stellar remnant, whose light curve can be used to measure 
cosmic distances. In archival Hubble Space Telescope imaging, three lensed images of the supernova 
are detected with relative time delays of <200 days. We predict a fourth image will appear close to the 
cluster core in the year 203742. Observation of the fourth image could provide a time delay precision 
of ~ 7 days, < 1% of the extraordinary 20 year baseline. The SN classification and the predicted 
reappearance time could be improved with further lens modelling and a comprehensive analysis of 


systematic uncertainties. 


We discovered AT 2016jka using data from the Hub- 
ble Space Telescope (HST) program REsolved QUIEs- 
cent Magnified Galazies (REQUIEM, HST-GO-15663, 
PI:Akhshik) [6] (see a summary of observations in Sup- 
plementary Table 1). The REQUIEM project targets 
massive galaxies with low specific star formation rates 
that have been magnified by strong gravitational lens- 
ing. The brightest and most spectacular galaxy targeted 
by REQUIEM is MRG-M0138, a massive red galaxy 
(MRG) at z = 1.95 [ref 7] behind the galaxy cluster 
MACS J0138.0-2155 [ref 8]. MRG-M0138 is quadruply 
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lensed by a foreground galaxy cluster at z = 0.338. Dur- 
ing analysis of observations obtained 13-14 July 2019 
we discovered three point sources that were present in 
archival HST images from 18-19 July 2016, part of the 
program (HST-GO-14496; PI:Newman) that first con- 
firmed the MRG-M0138 galaxy as a strongly-lensed ob- 
ject (Methods: Observations). Each point source is 
within 5 arcseconds of one of the four MRG-M0138 im- 
ages. None of the three point sources are present in 
the REQUIEM HST data in 2019 (Fig. 1). We infer 
that these are multiple images of a single astrophysical 
transient in MRG-M0138, most likely a SN. 

To construct a lens model for the MACS J0138.0-2155 
cluster we use the LENSTOOL software [9; 10] (Methods: 
Lens Modelling). To avoid unintended bias, we kept the 
lens model development completely separate from the 
analysis of the SN. Only upon completion of both were 
the results combined for the analysis described here. 
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Figure 1. Overview of the MACS J0138.0-2155 cluster field and the AT 2016jka discovery. The wide-field view in a) 
is 40” on a side with ticks indicating 10” intervals. Labels indicate the locations of the observed AT 2016jka images (SN1-SN3), 
the expected future images (SN4-SN5), and the multiply-imaged host galaxy (H1-H4). Labels 3.1-3.4 indicate the locations of 
a separate multiply-imaged [O11]-emitter at z = 0.77 used to help constrain the lensing potential. Panels b—i show 4” cutouts 
around the lensed SN images with 1” ticks. Panels bcde show the imaging from 2016 July where the SN was visible and fghi 
show the later imaging from 2019 July where the the SN has faded away. Panels aei indicate the location of the fourth image 
predicted to appear in +2037, with an ellipse demarcating the 68% probability region predicted from the LENSTOOL model. Panel 
a shows the expected location of the final and highly demagnified fifth image (“SN5”). The three-colour images are generated 
from the WFC3/IR filters with R=F105W or F110W, G=F125W and R=F160W or F140W (as indicated in panels dh). All 
panels use the late-epoch F125W imaging for the green channel; nevertheless, it is immediately clear that the SN2 image is 
substantially bluer than the other two, which helps to constrain the relative age of each SN image and the transient classification 


as a likely Type Ia supernova explosion. 


The input model constraints are the positions and red- 
shifts of the MRG-M0138 galaxy at z = 1.95 (both the 
galaxy’s centroid position and the SN location in each 
image) as well as a multiply-imaged background galaxy 
at z = 0.766, both having secure spectroscopic redshifts 
(Supplementary Table 3). We model the mass distribu- 
tion in the cluster core as the combination of a cluster- 
scale and galaxy-scale potentials (Supplementary Table 
4; Extended Data Figure 1). From this model we derive 
estimates for the lensing magnification and time delay 
of each of the SN images, including two predicted fu- 
ture images (Table 1). The lens model predicts that the 
SN should appear in the fourth MRG-M0138 image in 


the year 203742, demagnified with u = 0.4 + 0.2. A 
fifth image will also appear at a still later date, located 
near the center of the cluster and much more signifi- 
cantly demagnified, so it will not be easily observable. 
We anticipate that future lens modelling of the cluster 
will improve on these predictions primarily by exploring 
a wider range of mass models and incorporating more 
observational constraints (Supplementary Note: Future 
Work). For example, our analysis adopted only a single 
form for the density profiles, and did not incorporate 
constraints from stellar kinematics or pixel-level surface 
brightness data from the multiply-imaged systems. Al- 
though our LENSTOOL model does a good job of repro- 
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Table 1. SN lensing observables. 


Image R.A.* Dec.* Age? Ale Atprea® Hpred 
[hh:mm:ss] [dd:mm:ss] [days] [days] [days] 

SN1 01:38:03.63 —21:55:50.38 9213 — — 51 

SN2 01:38:02.96 —21:55:47.26 —24+16 1143$ 82 + 62 743 

SN3 01:38:02.42 —21:55:38.47 107436 —17+!13 -19434 3.90.5 

SN4 01:38:04.15+0.36 —21:55:24.73+0.43 — — 7742 + 5407 0.4+0.2 

SN5 01:38:03.77+0.1  —21:55:31.74+0.1 — = 9463 +770° < 0.01 


a Coordinates are given in the J2000 reference frame, as measured for images 1-3 and as predicted for 
images 4 and 5. Uncertainties in the predicted position of images 4 and 5 are in arcseconds. 

P The measured age of each SN image is the number of days relative to peak brightness, in observer-frame 
days. A negative value means the image was observed prior to peak brightness. 


© Time delays are reported relative to image SN1. 


d The combination of the models for the light curve evolution and the lens time delays implies a date of 


2037 + 2 for the peak brightness of Image 4 and 20424 


ducing the morphology of the host galaxy images H1-H3, 
it does not reproduce image H4 as well (Supplementary 
Note: Host Image Morphology Comparison). 

If we can estimate the age of each SN image then we 
can derive direct observational constraints on the rel- 
ative lensing time delays. For this goal, it is helpful 
to have a firm determination of the transient’s class. 
Expected time delays and magnifications from the lens 
model exclude any of the various rapidly evolving and 
low-luminosity stellar transient classes, strongly suggest- 
ing that it is a SN. The first-order SN distinction re- 
maining is between a Type Ia SN—the explosion of a 
white dwarf star in a binary system—and a core col- 
lapse SN (CCSN)—the end-point of a star with mass 
> 10Mọo. The properties of the host galaxy can inform 
this classification because CCSN are limited to galaxies 
with young stellar populations. Limits on the specific 
star formation rate and age for this host, MRG-M0138, 
show it to be a massive but very quiescent and evolved 
galaxy, unlikely to retain any significant population of 
high-mass stars |7]. Based on observed properties of the 
host galaxy alone, we find a 62-75% probability that 
AT 2016jka is a Type Ia SN (Methods: Classification). 

Adopting the lens model magnifications for the three 
observed SN images (Table 1) we can locate each SN 
image in colour-magnitude space (Fig. 2). After mag- 
nification correction, all three images are still brighter 
than expected for a Type Ia SN, which may indicate that 
a lens-modelling degeneracy is at play. Nevertheless, 
the magnification-corrected AT 2016jka data are more 
consistent with the Type Ia population than any CC 
SN sub-class (Fig. 2a,b and Extended Data Fig. 2). 
AT 2016jka also demonstrates the expected evolution of 
a Type Ia SN colour and brightness over ~ 100 days 
(Fig. 2c). By also including the model-predicted time 


t4 for Image 5 (which will likely be undetectable). 


delays, we can treat the three SN images effectively as 
three points on a common SN light curve, and we find 
p(Ia)=94% (Methods: Classification). This composite 
light curve is shown in Figure 3, with the best-fitting 
Type Ia SN model. Extended Data Fig. 3 also shows a 
random draw of light curves from the Monte Carlo sam- 
pling for all three major SN sub-classes. For the remain- 
der of this analysis, we proceed under the assumption 
that AT 2016jka is indeed a Type Ia SN. An improved 
classification could be achieved with spectroscopy and 
multi-band photometry upon arrival of the fourth im- 
age. In that case the analysis that follows here could be 
revised to achieve similar results with a different under- 
lying SN model. 

The colour of a Type Ia SN evolves substantially over 
its lifetime as the photosphere expands and cools, reveal- 
ing different layers of the expanding shell and driving 
episodes of recombination [11]. Since the phenomenon of 
gravitational lensing in general is achromatic, this colour 
evolution makes it possible to derive an age constraint 
that is largely independent of the lens model (Methods: 
Colour Curve Age Constraints; Extended Data Fig. 4 
and 5). Combining this information with magnification 
constraints from the lens modelling helps break param- 
eter degeneracies, yielding measured delays in Table 1 
(Methods: Light Curve + Lens Model Age Constraints; 
Extended Data Fig. 6 and 7). 

Using these measured time delays, we created a re- 
constructed form of the intrinsic light curve and colour 
curve of AT 2016jka, shown in Fig. 3. Remarkably, the 
ages of image 1 and 3 are constrained to better than +20 
days, despite having only a single epoch of photometric 
data. These uncertainties may be further reduced when 
the future fourth image is observed with high-precision, 
multi-epoch photometry. Such a light curve will pin 
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Figure 2. Classification information for AT 2016jka based on its position in colour-magnitude space. (a) 
The observed photometry for the three SN images is shown as open markers. Vertical dotted lines show the lens magnification 
corrections using our preferred lens model, with a filled grey marker indicating the corrected (demagnified) magnitude. Error bars 
indicate the observational and systematic uncertainty, including the range of alternative magnification corrections encompassed 
by lens model variations. Contours show the population distributions for normal SNe of Type Ia (red), Type Ib/Ic (gold), and 
Type II (green), drawn with two contour levels enclosing 68% and 95% of each SN population. Each SN sub-class was simulated 
at z = 1.95, and samples from their expected light curves were drawn uniformly in time. (b) Marginalized distributions along 
the colour dimension for the three SN sub-classes (using the same colour scheme). The simulated populations have been scaled 
according to the expected explosion rates in the SN host galaxy, based on its stellar population properties, so the y-axis scaling 
is effectively a relative probability. Open markers show the observed colours again, at arbitrary y position. (c) Zoomed-in view 
of the colour-mag space marked by the grey box in panel a. The evolution of a typical Type Ia SN at z = 1.95 is shown by 
a coloured line, with the line colour indicating SN age in observer-frame days relative to peak brightness. White diamonds 
correspond to the times labeled on the colourbar below. Grey shading shows the typical range of luminosities and colours 
observed for the Type Ia SN population in the nearby universe. Although the magnification-corrected data are brighter than 
expected for most Type Ia SN, they are consistent both with the overall Type Ia SN population, and with the Type Ia SN 
colour-mag vs time curve. 


down the intrinsic SN light curve parameters that are 
shared by all images, and break remaining parameter de- 
generacies. Improvements to the lens modelling will also 
be essential, to better estimate and minimize systematic 
biases that may arise from the necessary magnification 
and time delay corrections. 

AT 2016jka and other lensed SN like it could eventu- 
ally contribute to mapping the cosmic expansion his- 
tory and measuring the effects of dark energy, which 
appears to be driving an accelerating cosmic expansion 
rate [12; 13]. Recent investigations into the expansion 
rate of the universe (the Hubble-LeMaitre constant; Ho) 
have found that measurements from the local universe 
are significantly different from the value inferred from 
measurements of the cosmic microwave background ra- 
diation [14; 15]. The community is actively attempt- 
ing to resolve this “Ho crisis” by mitigating systematic 


uncertainties or discovering new physics from the early 
universe [16]. Either resolution will require multiple in- 
dependent cosmological probes. In recent years, lensed 
quasar time delays have provided a valuable independent 
tool for this, with seven high-precision measurements to 
date [17]. As the sample of time delay lenses grows to 
several dozen (including lensed SNe like AT 2016jka), it 
is expected to deliver a measurement of Ho with 1% 
precision [18]. 

Looking beyond the Ho crisis, determining the na- 
ture of dark energy and how it may evolve over time 
is a primary goal for the large-scale cosmology experi- 
ments of the 2020s [19; 20]. A future sample with ~ 100 
well-measured lensing time delays would be a compet- 
itive tool for dark energy studies [4; 18]. Events like 
AT 2016jka could be an important part of this time- 
delay cosmology sample, but to date there have been 
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Figure 3. The reconstructed light curve and colour curve for AT 2016jka. Panel (a) shows the light curve and panel 
(b) shows the colour curve, both incorporating lens model corrections for the time delay and magnification. The grey shaded 
region covers the 68% confidence interval of the best-fit Type Ia model (using the SALT2 model [31]), with the median model 
shown as a solid line. Observed photometric data are shown as coloured markers. The x-dimension error bars on each data 
point represent the lens model time delay uncertainties (Table 1). The y-dimension uncertainties in both panels incorporate 
the photometric uncertainty. For the light curve (a), the y error bars also include uncertainty in the lens model magnification 
(Table 1). The colour measurements (y values, panel b) do not require lens model correction. Parameters for the best-fit SALT2 
models in grey are shown in the legend, and were obtained using the joint posterior of the colour curve and light curve methods 
(described in Methods:Time Delay Estimation; shown in Extended Data Fig. 7). 


only two lensed SN observed with multiple images. The 
first, SN Refsdal, was a peculiar Type II SN whose im- 
age with the longest delay was missed [21; 22]. The 
second, SN 2016geu, was a Type Ia SN with short de- 
lays that make high precision time delay measurements 
impossible [23; 24]. An earlier discovery of an unusu- 
ally luminous SN was also shown to be a strongly-lensed 
Type Ia SN [25], though the multiple images were not re- 
solved. This makes AT 2016jka just the third discovery 
of a lensed SN resolved into multiple images. 

Future large-scale surveys such as the Vera C. Rubin 
Observatory and the Nancy Grace Roman Space Tele- 
scope will observe dozens to hundreds of lensed SN over 
their mission lifetimes. The vast majority of these will 
be lensed by galaxy-scale deflectors [26; 27] and thus 
will have significantly shorter delays (of order 10-100 
days). Since it is the fractional time delay uncertainty 
that propagates through to any time delay distance mea- 
surement, the extraordinarily long time delays of cluster- 
lensed SNe like AT 2016jka can deliver significantly bet- 
ter time delay precision, with comparable observational 
cost (Supplementary Note: Future Discoveries). In fact, 
the long time baseline of lensed SN like AT 2016jka effec- 


tively insures that their cosmological precision is not lim- 
ited by time delay measurement uncertainty. Cluster- 
scale lenses are more complex than galaxy-scale lenses, 
but they generally have “independent” measurements of 
the magnification from several multiply-imaged systems 
in the same field. Modelling cluster lenses is very differ- 
ent from galaxy lenses, so objects like AT 2016jka can 
provide a valuable check on systematics for the larger 
sample of transients used in time-delay cosmology. 

The first multiply-imaged SN discovery, SN Refsdal, 
has shown that time delay cosmography with a cluster- 
lensed SN is viable [28; 22]. However, the observing 
campaign for SN Refsdal was extraordinary, deploying 
more than 75 HST orbits over three years. Significant 
observational investment has also been required for high- 
precision time delay measurement of lensed quasars, 
such as decade-long programs [29] or daily monitor- 
ing for high-cadence light curves [30]. In the case of 
AT 2016jka it will be possible to achieve similar time 
delay precision over the 20-year baseline with just a sin- 
gle imaging epoch as the anchor point. A sample of 
AT 2016jka-like events could be developed with regular 
monitoring of cluster-scale lenses, partnered with mod- 
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est follow-up to characterize any lensed SN discovered 
(Supplementary Note: Future Work). 

HST observations enabled us to find this SN. We an- 
ticipate that HST may be de-orbited and make its final 
plummet to Earth around the time of the reappearance 
of AT 2016jka, so we coin the name SN Requiem as an 
ode to the vast new discovery space that HST continues 
to unveil. 
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Extended Data Figure 1. Elements of the MRGO138 cluster lens model. The model comprises 37 potentials in 
total: the BCG (red), 32 cluster members (yellow), three perturbers (cyan), and the main cluster potential (pink). Labeled x 
symbols indicate the positions of the SN, host, and one additional multiply-imaged galaxy with a secure redshift used as model 
constraints (Supplemental Table 3). The filters used to generate the color image are as in Fig. 1, and tick marks are separated 


by 10 arcsec. 
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Extended Data Figure 2. The position of AT 2016jka in color-magnitude space. Colored points show simulated 
photometry for normal SNe of Type Ia (red), TypeIb/Ic (gold), and Type II (green), with 10,000 simulated SN in each sub-class 
(not all apparent on this plot). Histograms above and below show the marginalized distributions that have been rescaled to 
represent posterior probability density functions. They are normalized to integrate to unity, then multiplied by the SN sub- 
class priors based on the host galaxy stellar population (row b in Supplementary Table 2). Open markers show the observed 
photometry of the SN. Dotted vertical lines mark the magnification correction based on the preferred LENSTOOL model 
(model E, described in Methods: Lens Modeling). Closed markers show the resulting magnification-corrected photometry, with 
asymmetric error bars reflecting the systematic uncertainty derived from the five lens model variants. Horizontal error bars 
in the upper panel indicate the observed uncertainty in the SN color (not affected by lensing). The relevant SN photometry 
markers are repeated in the histogram side-panels with arbitrary vertical positions. All three SN images are located in regions 
of color-magnitude space that are expected to be dominated by Type Ia SN. 
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Extended Data Figure 3. A representative set of light curve and color curve models from the STARDUST2 
classification algorithm. Panels a-f show F105W and F160W, as indicated, plotting the model light curves in black and 
photometry as red markers. Panels g-i show the F105W-F160W color curves and color data. All data points as shown have 
been corrected for magnification and shifted in time using the preferred LENSTOOL model (model E, described in Methods: 
Lens Modeling). Plotted error bars include the measurement uncertainty and the lens modeling magnification uncertainty. Data 
points in the right column also include this magnification uncertainty, even though cluster-scale lensing is achromatic, because 
the STARDUST?2 analysis was done on the light curve data, not the color data directly. In all panels the first data point is SN 
image 2, followed by image 1 and image 3. In each panel the black curves show 200 SN light curve models drawn at random 
from the nested sampling sequence of the STARDUST2 (sncosmo) classification. 
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Extended Data Figure 4. Color-based age constraints for AT 2016jka. Constraints are shown separately for Image 
SN1 (panels a and b), image 2 (c and d), and image 3 (e and f), using the methodology described in Methods: Color Curve Age 
Constraints. Large lower panels (b, d and f) show the observations and model fits. Each magenta shaded region shows the lo 
range of the measured F105W-F160W color, which corresponds to a U-V color in the rest-frame. The model fits are shown as 
grey shaded regions, indicating the 68% confidence interval of the best-fit SALT2 color curve, with the median model shown as 
a solid line. The small upper panels (a, c and e) show the posterior for the age of each image from SNTD, using a prior on the 
SALT2 color parameter (c) based on known population characteristics of SNla. The effect of adding this prior is slight, with no 


significant deviation from the best-fit value of c = 0.02*9'}3). 
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Extended Data Figure 5. Marginalized and joint posterior distributions for the color curve age constraints 
measured in this analysis. Two-dimensional plots show MCMC sampling points as discrete dots, with contours for the high 


density regions, drawn at 0.5, 1, 1.5 and 
dashed vertical lines at the mean and 4 


20 levels. Marginalized (1D) distributions are shown at the top of each column, with 
tlo (marking 16%, 50% and 84% levels). We use a weak prior on the SALT2 color 


parameter (c), and set the SALT2 stretch parameter (x1) to 0. This method is fully independent of lens modeling. The table 
in the upper right lists all priors, observations, and lens model information used for SN age estimates in this work. Only the 


highlighted components were used for th 


e constraints shown here. 
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Extended Data Figure 6. Light-curve-based age constraints for AT 2016jka. Constraints are shown separately for 
Image SN1 (panels a and b), image 2 (c and d), and image 3 (e and f), using the methodology described in Methods: Light 
Curve Age Constraints. Large lower panels (b, d and f) show the observations and model fits. Each orange shaded region shows 
the lo range of the measured F160W magnitude after lens model correction (see Table 1), which corresponds to roughly V band 
in the rest-frame. The model fits are shown as grey shaded regions, indicating the 68% confidence interval of the best-fit SALT2 
light curve, with the median model shown as a solid line. Small upper panels (a, c and e) show the posterior distributions from 
SNTD for the age of each image that is independent of the lens model (magenta, same as Extended Data Figure 7), using the 
preferred lens model E (light blue), and the combination of both methods (orange). 
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Extended Data Figure 7. Marginalized and joint posterior distributions for the final age constraints measured 
in this analysis. Two-dimensional plots show MCMC sampling points as discrete dots, with contours for the high density 
regions, drawn at 0.5, 1, 1.5 and 20 levels. Marginalized (1D) distributions are shown at the top of each column, with dashed 
vertical lines at the mean and +1o (marking 16%, 50% and 84% levels). We use the color curve posterior as the prior for 
light curve fitting with lens model E, and include weak priors on the absolute magnitude of a SNIa (MB) and the SALT2 color 
parameter (c), and set the SALT2 stretch parameter (21) to 0. 
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Methods 


OBSERVATIONS 


The observations used in this work are summarized 
in Supplementary Table 1. We processed all HST ob- 
servations using the Drizzlepac software utilities [32], 
aligned to a common astrometric reference frame and re- 
sampled to 0.1 arcseconds per pixel. We then identified 
isolated and unsaturated stars in each image and used 
them to create an effective point spread function (ePSF) 
with 4x oversampling, using the photutils package 
from the astropy software suite [33]. To measure the 
SN photometry we followed two tracks. As our pri- 
mary method we performed ePSF fitting directly on the 
F105W (Y band) and F160W (H band) images where 
the SN was apparent. This fitting allowed for a constant 
background flux to account for both the sky brightness 
and the background light of the cluster and host galaxy. 
As a second approach, we created “pseudo-difference im- 
ages” by re-scaling the F110W and F140W images col- 
lected in 2019 (in which the SN is not present). The 
transmission functions of the F110W and F140W fil- 
ters are broader than F105W and F160W, and do not 
strictly overlap in wavelength. The optimal scaling fac- 
tor to produce a clean subtraction therefore depends 
on the spectral energy distribution of the source. We 
set the scaling to 0.62 and 1.17 for F110W-to-F105W 
and F140W-to-F160W, respectively. These values pro- 
duced visually clean subtractions of the SN host galaxy 
MRGO0138—meaning that they minimize the residual 
flux in the pseudo-difference images. We then performed 
ePSF fitting on the SN in each pseudo-difference image, 
as before. Both sets of photometry agree to within one 
standard deviation. The reported photometry in Sup- 
plementary Table 2 are the measurements from the first 
method (collected directly from the un-subtracted im- 
ages). Also reported in Supplementary Table 2 are flux 
densities and uncertainties measured within D = 0.”7 
circular apertures at the position of the SN3 image in 
the 2019 HST visits where the SN is not detected. 


VLT SPECTROSCOPY 


We make use of integral field spectroscopic data 
obtained on the cluster core of MRG0138 with 
VLT/MUSE, publicly available as part of the program 
0103.A-0777(A) (PI: Edge). Three exposures of 970 sec 
each were taken with a small dithering offset and 90 
degree rotations in between. This dataset was reduced 
and analysed using the MUSE data reduction pipeline 
v.2.7 [ref 34] for basic calibration (bias, flat-field, wave- 
length, LSF, geometry) as well as flux calibration, sky 


subtraction and astrometry. We also make use of the 
self-calibration technique [35] to remove illumination 
systematics, specifically tuned for the case of crowded 
fields in the central region of galaxy clusters (Richard 
et al. in prep.). The combined datacube is then pro- 
cessed through ZAP [36] which applies a PCA technique 
to remove sky subtraction residuals. The final datacube 
covers the central 1x1 arcmin? around the cluster center 
with 0.2”x0.2”1.25A pixels. 

We have extracted spectra for each HST detected 
source and inspected them for redshift measurements. 
In addition, we have run the muselet software (part 
of the MPDAF package [37]; mpdaf.readthedocs.io) to 
search for line emitters not directly associated with HST 
sources [38; 39]. Apart from the lensed quiescent galaxy, 
we measured spectroscopic redshifts for cluster members 
and one ring-like background galaxy north of the BCG 
at z = 0.766. Finally, we have measured the velocity 
dispersion of the BCG to be 390410 km s~!. This mea- 
surement was not available for use in our blind lens mod- 
elling, but could provide a useful constraint for future 
lens model improvements (Supplementary Note: Future 
Work). 


LENS MODELLING 


To model the mass distribution in the MACS J0138.0- 
2155 cluster core we use the latest version of LENSTOOL 
[9] (git-cral.univ-lyon1.fr/lenstool), which performs a 
Bayesian analysis with an MCMC sampler to estimate 
the best fit and uncertainty on each parameter of the 
mass distribution. 

The strong-lensing constraints used are the locations 
of multiple images found in HST and MUSE/VLT. We 
group them into 3 systems: (a) the 4 images of the 
quiescent galaxy hosting MRG0138-SN, (b) the 3 ob- 
served images of AT 2016jka assumed to be at the same 
redshift, and (c) the diffuse arc-like structure identified 
in HST and confirmed as an [Ou]\3727 emitter in the 
VLT/MUSE datacube (see previous section). The im- 
age coordinates and redshifts used in the lens model are 
summarised in Supplementary Table 3. 

The cluster mass modelling is performed similarly to 
other massive strong lensing clusters observed with HST 
[40]. The total mass distribution is parametrized as a 
combination of multiple dPIE (double Pseudo Isother- 
mal Elliptical) profiles describing both cluster-scale and 
galaxy-scale dark matter haloes. dPIE are elliptical 
isothermal profiles with both a core radius and a cut ra- 
dius where the density flattens and drops, respectively. 
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In the case of MRGO0138 the mass distribution is dom- 
inated by a single mass concentration centered on its 
Brightest Cluster Galaxy (BCG). We therefore use a 
single cluster-scale halo at a fixed cut radius of 1 Mpc. 
We add a single galaxy-scale halo on each cluster mem- 
ber, where the shape parameters (halo center, ellipticity) 
are fixed to their measured HST morphology and their 
core radius is negligible (fixed at 0.15kpc). Extended 
Data Fig. 1 shows the locations of all components of the 
cluster model, including 32 cluster members. Cluster 
members were identified by the combination of red se- 
quence selection (based on the F814W—F160W colour) 
and MUSE spectroscopy. 

The majority of cluster members are elliptical galax- 
ies selected from the red sequence, and to reduce the 


number of parameters we assume they follow the scaling 
1, \a/4) 
L* 


relations: o = o* for the velocity dispersion, 


(1/2) 
and Tout = Tout (4) , assuming the Faber-Jackson 


relation and a constant M/L ratio respectively. o* and 
re, are model parameters for a cluster member at the 
characteristic luminosity L*. Following the discussion 
in [41] we fix o* = 158 km/s and r¢,,=45 kpc. 

We individually optimise the ø and rey; parameters 
for 4 specific galaxies which are not expected to fol- 
low the aforementioned scaling relations: the BCG and 
three perturbers P1 to P3. These perturbers are ei- 
ther blue gas-stripped galaxies infalling into the clus- 
ter core, and/or located very close to the images of the 
AT 2016jka host galaxy, perturbing its apparent mor- 
phology with additional lensing. This choice of per- 
turbers is similar to the ones used in the model by [7]. 

LENSTOOL optimises the parameters of the model by 
minimising the overall root mean square dispersion 
(RMS) between the predicted and observed locations of 
the multiple images. The best fit parameters of each 
mass component are provided in Supplementary Ta- 
ble 4. Uncertainties reported there are derived from the 
MCMC models, by sampling their posterior probability 
distribution. 

We developed five lens model variants blindly (i.e., 
without knowing the impact of each lens model varia- 
tion on the transient classification or time delay infer- 
ences). Model A was the first viable model developed, 
which did not include additional perturbers, and did 
not include the additional lensed background source at 
z = 0.7663. In Model B we allowed for the location 
of the main cluster dark matter halo to be free, with 
an offset from the reference position taken at the BCG 
center. Model C allowed the same central position off- 
set and also relaxed the constraints on the BCG (o and 
Teut)- Model D fixed the primary dark matter halo at 


the BCG center, but still relaxed the constraints on the 
BCG o and reut. The final model, and the one selected 
as the preferred model prior to unblinding, is model E, 
which includes all four perturbers described above, and 
includes the additional background object at z = 0.7663. 
Note that all of the lens model variants A-D would give 
a slightly larger magnification for all three SN images. 
This means that our estimated systematic uncertainties 
are one-sided, as can be seen on Fig. 2 and Extended 
Data Fig. 2. 

This model is then used to predict the magnifica- 
tion and time delays for the three observed images of 
AT 2016jka, as well as the location of the fourth and fifth 
images, which are still to appear. These predictions are 
summarised in Table 1, but note that the magnification 
predictions for image SN4, as well as all predictions for 
SN5, should be treated with caution, as the lens models 
presented here have known shortcomings. For example, 
the best-fit velocity dispersion for the BCG is 700 km/s 
(Supplementary Table 4), though we have measured this 
property from MUSE spectra to be 390+10 km/s. Fur- 
ther lens modelling is needed to incorporate such ad- 
ditional constraints and to fully quantify potential sys- 
tematic biases (Supplementary Note: Future Lens Mod- 
elling). 

The final LENSTOOL model E reproduces all multiple 
system positions with an RMS of 0.15”. The lo un- 
certainty from model E for the SN4 location is the el- 
lipse overlayed on panels e and i in Fig. 1. As the 
lens model reproduces the location of the SN images 
within a small uncertainty, these predictions are com- 
puted with LENSTOOL using the barycenter of all source 
positions corresponding to images SN1, SN2 and SN3 
as the same reference source position. Comparing our 
lens modelling to the previous model of this cluster from 
[7], we find that the magnification estimates are broadly 
consistent, though systematically lower (Supplementary 
Note: Comparison to Previous Lens Modelling). 


CLASSIFICATION 


The lens modelled time delays between the images are 
~ 100 observer-frame days, but we see that three images 
of the transient are visible simultaneously. From this we 
can infer that the visibility time of the transient in the 
z = 1.95 rest-frame must be at least ~30 days. Sim- 
ilarly, with expected magnifications in the vicinity of 
u ~ 10, the measured apparent magnitudes near 23 AB 
mag translate to a rest-frame absolute magnitude near 
Mpg ~ —19.5 mag (too bright to be a nova, luminous 
blue variable, or other low-luminosity stellar transient). 
Taken together, these indicators strongly suggest that 
the transient is a supernova (SN). 
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Although we have invoked the lens model in this anal- 
ysis, we note that the inferences are not strongly depen- 
dent on the specific lens model predictions. To make 
the observed transient images consistent with a fast or 
low-luminosity transient, the time delays and/or mag- 
nifications would have to be changed by more than a 
factor of 2. In the analysis to follow, we will work under 
the assumption that AT 2016jka is a SN. 


SN SUB-CLASSIFICATION BASED ON HOST GALAXY 


With this transient identified as a SN, we now seek 
to identify the most likely SN type, under the assump- 
tion that it belongs to one of the three most common 
sub-classes (Ia, II, Ib/c). We first use two methods that 
rely only on measured properties of the host galaxy to 
circumstantially infer the type. This inference is less 
strongly dependent on the lens model, helping to reduce 
any bias associated with a lens model-dependent classi- 
fication. 

Although Type Ia SNe are found in all types of galax- 
ies, CCSNe are limited to galaxies with relatively young 
stellar populations. We can therefore infer some infor- 
mation about the SN type using the observed host prop- 
erties combined with knowledge of the relative rates of 
Type Ia and CCSNe in different stellar populations [42]. 
In the case of the host galaxy MRGO0138, we have a 
very well-constrained spectral energy distribution (SED) 
extending out to far-infrared wavelengths with Spitzer 
IRAC data [7; 43]. From the SED fitting we derived the 
host galaxy’s rest-frame B-— K colour and absolute mag- 
nitude, Mx, which serve as proxies for the stellar pop- 
ulation age and have been empirically calibrated with 
SN rates in the local universe [44]. We adopt a lensing 
magnification correction using LENSTOOL model E to get 
Mx. The B—K colour is not affected by the foreground 
lens. Using the galsnid method [44] we derive a 75% 
probability that AT 2016jka is of Type Ia (Supplemen- 
tary Table 5). We used host galaxy image 2 for this 
purpose, with the flux-weighted harmonic mean magni- 
fication u = 8.3 to derive Mx (see Supplementary Ta- 
ble 6). We also evaluated the other host galaxy images 
and found no change in the resulting SN classification 
probability. 

As an alternative host galaxy classification constraint, 
we use the derived properties of the host galaxy stel- 
lar population directly, rather than adopting colour and 
magnitude proxies. The MRG-M0138 galaxy has high 
mass (logi9(M/M.) = 11.7) but is a very quiescent 
galaxy, with a specific star formation rate of ~ 107113 
yr! and a stellar population that is well-matched by 
an exponential star formation history with an age of 1.4 
Gyr [ref 7]. The massive stars that end as CCSN explo- 


sions have main-sequence lifetimes of < 40 Myr [ref 45], 
making it unlikely that CCSN progenitors make up a sig- 
nificant fraction of the MRG-M01388 stellar population— 
though the high total stellar mass makes it possible that 
pockets of young stars are present. We define classifi- 
cation probabilities based on the projected SN rate for 
each SN sub-class, derived from the host galaxy’s stellar 
mass and star formation rate [46]. This yields a 62% 
probability that AT 2016jka is of Type Ia (Supplemen- 
tary Table 5). 


SN SUB-CLASSIFICATION BASED ON SN 
PHOTOMETRY 


To improve the classification of the AT 2016jka sub- 
type, we now bring in observed photometry of the SN it- 
self, and again we adopt two methods. The first method 
uses only magnification information from the lens model, 
and the second uses both the modelled magnification 
and time delay predictions. In both cases we adopt the 
stellar-population-based host galaxy classification prob- 
abilities as priors. 

Extended Data Fig. 2 illustrates the first approach. 
After applying the magnification corrections, each of the 
three images of the SN are mapped to colour-magnitude 
space. We then treat each observed point (each SN im- 
age) separately, comparing their colour-magnitude loca- 
tion to a simulated population of unlensed SNe. Our 
simulation uses the sncosmo package [47] to generate 
10,000 SN for each of the three principal SN sub-classes 
(Ia, Ib/c, II), all at z=1.95. We then compute the num- 
ber of simulated SN within a rectangular region around 
each observed point. The width of this sampling region 
is set to 3 times the observed colour uncertainty, and 
the height is equal to the lens-modelling magnification 
uncertainty. The u uncertainty used here includes an 
estimate of the systematic uncertainty for each SN im- 
age, derived from the spread of magnifications across 
lens model variants (similar to the methodology of [7]). 
Note, however, that the lens model variations evaluated 
here do not vary the assumption of the density pro- 
file, which can strongly affect magnifications (see also 
the Supplementary Note: Comparison to Previous Lens 
Modelling). 

We take the number of simulated SN for each type 
within the sampling region as an estimate of the likeli- 
hood that AT 2016jka belongs to that class. Note that 
in this case there is no need to apply a cut to the sim- 
ulated sample to account for detectability, because the 
5o limiting magnitude of our HST observations is 26.5 
AB mag, and after accounting for magnification of ~1.5 
mag (Table 1) this becomes mijn, ~ 28 mag. This 
means that all the points shown in Extended Data Fig. 2 
(and therefore all simulated SN entering our classifica- 
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tion counts) would be easily detectable in our HST imag- 
ing. Multiplying by the prior probabilities derived from 
host galaxy properties, we finally derive the probability 
that the SN is of Type Ia as p(Ia) = 0.92, 0.98, and 0.95 
from the three SN images SN1, SN2 and SN3, respec- 
tively. Supplementary Table 5 reports the mean of our 
three classification probabilities for each sub-class. 

As a second photometric classification of this SN, we 
used the STARDUST2 Bayesian light curve classification 
tool [48], which is also built on the underlying sncosmo 
framework. Here we adopt both the predicted mag- 
nifications and time delays from the best lens model, 
which allows us to put the photometry from the three 
images together as a composite “light curve” and com- 
pare against simulated light curves. STARDUST2 uses 
the SALT2-extended model to represent Type Ia SN 
[31; 49] and a collection of 42 spectrophotometric time 
series templates to represent CCSN (27 Type II and 
15 Type Ib/c). These CCSN templates comprise all 
of the templates developed for the Supernova Analy- 
sis software SNANA [50], derived from the SN samples of 
the Sloan Digital Sky Survey [51; 52; 53], Supernova 
Legacy Survey [54], and Carnegie Supernova Project 
[55; 56; 57]. With STARDUST2 we use a nested sampling 
algorithm to measure likelihoods over the SN simula- 
tion parameter space. Extended Data Fig. 3 shows the 
magnification- and time-delay-corrected photometry of 
AT 2016jka and a random sampling of light curve mod- 
els from the sncosmo nested sampling algorithm em- 
ployed by STARDUST2. Nested sampling is a Monte Carlo 
method that traverses the likelihood space in a manner 
that samples the Bayesian likelihood [58]. These sample 
light curves and colour curves therefore give a visual rep- 
resentation of how well each SN sub-class (Ia, II, Ib/c) 
can match the observed data. This figure shows that the 
limited photometric data can be reasonably well fit by at 
least one model from any of these three sub-classes. The 
density of curves in the top panels demonstrates that the 
Type Ia SN model is consistently a good match to the 
data. However, the range of model parameters that al- 
low such a fit to the data is much more limited for the 
heterogeneous CCSN types. To compute the posterior 
probability distribution we adopt priors for each of the 
three SN classes, again using the classification probabili- 
ties derived from the AT 2016jka host galaxy stellar pop- 
ulation properties (Supplementary Table 5). Marginal- 
izing the posterior probability distributions over all free 
parameters, we find a 94% probability that AT 2016jka 
is of Type Ia (Supplementary Table 5). 

The combination of evidence from the derived host 
galaxy properties and SN photometry supports the con- 
clusion that AT 2016jka is a Type Ia SN with > 90% con- 


fidence. Although we have adopted some lens modelling 
corrections for all of these methods, this conclusion is 
not sensitive to the choices we can reasonably make for 
modelling the lens. As shown in Extended Data Fig- 
ure 4, every LENSTOOL model we have evaluated locates 
the AT 2016jka de-magnified position within the region 
of colour-magnitude space dominated by Type Ia SN. 
Dust is also not a confounding factor here. The sim- 
ulations used in both SN-based classification methods 
include dust extinction at the source plane. If a signif- 
icant screen of dust exists in the lens plane, this would 
have the effect of making the SN appear dimmer and 
more red, so correcting for that extinction would move 
the AT 2016jka points upward and leftward on Extended 
Data Fig. 2, which would not shift it into the regions oc- 
cupied by Type II and Ib/c SN. 


TIME DELAY ESTIMATION 


We constrain the relative time delays of AT 2016jka 
by using two separate methods to estimate the age of 
the SN at each image during the single observed epoch. 
The preferred method of SN time delay measurements 
involves measuring the time of peak brightness for the 
SN at each image by fitting the light curves, and taking 
the difference between each measurement as the rela- 
tive time delay [59; 24; 60]. With only a single observed 
epoch, this method is impossible due to model parame- 
ter degeneracies, and we must rely on colour and bright- 
ness to constrain the age of each image of the SN. Such 
age estimates are sometimes referenced to the time of 
explosion, but in this case we use the observer’s conven- 
tion, setting age=0 as the time of peak brightness in the 
rest-frame B band (A ~ 4500 A). Each of these images 
stems from the same SN explosion, so the difference be- 
tween the measured age of each image is also a measure 
of the relative time delay. 

In all of the light curve fitting exercises described 
below, we also fix the SALT2 “stretch” parameter at 
x, = 0. This parameter defines the shape of the SN 
Ia light curve (the rate of decline in brightness). If we 
allow x; to be a free parameter, there is no useful con- 
straint on it. It is highly degenerate with the time delay 
between the images, which is of course a free parame- 
ter in all the fits. As a check, we have also tried fixing 
zı to other values from —1 to +1, and the time delay 
results change by less than 5 days, which is well within 
all error bars. Fixing x, in this way is comparable to 
the analysis that will be possible in the 2030s when the 
fourth SN image is observed. We expect that a SALT2 
fit to a well-sampled light curve from the fourth image 
will provide a tight constraint on xı. That measured zı 
will then be propagated back as a fixed parameter (with 
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small uncertainty) into revised fitting of the SN images 
1-3. Because of this, we do not incorporate the +5 day 
systematic uncertainty in the time-delay errors reported 
in Table 1. 


COLOUR CURVE AGE CONSTRAINTS 


We first attempt to constrain the relative time delays 
using the colour of each observed image, which is inde- 
pendent of the lens model and possible because the phe- 
nomenon of gravitational lensing is intrinsically achro- 
matic. One important caveat to this principle is that 
microlensing effects are not generally achromatic, be- 
cause the microlensing caustics may cause differential 
magnification on the scale of the SN radius [61; 62; 63]. 
Hence, if the expanding SN shell has a colour gradi- 
ent then microlensing may introduce spurious features 
in the observed colours of the SN [64; 65]. Ref. [61] 
found that such chromatic microlensing is most likely 
not present for lensed Type Ia SNe in the period up to 
about 25 rest-frame days after explosion (~15 observer 
days after peak brightness for AT 2016jka). Only im- 
age 2 is likely in the achromatic microlensing phase, but 
Ref. [61] found extremely small deviations in the rest- 
frame U — V colour curve due to microlensing at the 
68% confidence interval, and up to a ~ 0.2 — 0.4 mag 
difference with 99% confidence. While such extreme mi- 
crolensing could alter the results for images 1 and 3, it 
would not alter the measurement of image 2 as it is likely 
in the achromatic phase. Fortuitously, images 1 and 3 
are minima, which are less susceptible to high deviations 
owing to microlensing when compared to image 2, which 
is a saddle. 

We use version 2 of the SuperNova Time 
Delays (SNTD) package (publicly available at 
github.com/jpierell4/sntd) with documentation at 
sntd.readthedocs.io, which has several improvements 
over the original SNTD package [59]. The SNTD pack- 
age employs a nested sampling algorithm within three 
separate methods to measure time delays, and is de- 
signed to fully utilize the information present in SN light 
curve templates [66; 31; 67; 49] to reduce the impacts 
of microlensing and make more accurate measurements. 
We use the “colour” method present in SNTD, which 
attempts to reconstruct the intrinsic colour curve us- 
ing the SALT2 model as a template [31]. This method 
fits the age of each image simultaneously, while also 
varying the SN model parameters. This means we are 
finding a single set of SN model parameters to describe 
the intrinsic photometric evolution of the SN, and also 
finding the age (time from peak brightness) for each of 
the three images. 


The colour curve constraints resultiing from this pro- 
cess are shown in Extended Data Fig. 4. Joint and 
marginalized posterior distributions from SNTD for the 
SALT2 SN model parameters and the measured ages for 
each image are shown in Extended Data Fig. 5. In Ex- 
tended Data Fig. 4, one can see the measured colours 
intersect the model at two distinct locations for images 
2 and 3 of AT 2016jka—meaning there are two plausible 
ages. This results in a double peaked posterior distri- 
bution in Extended Data Fig. 5. This is caused by a 
model parameter degeneracy that could be broken in a 
way independent of the lens model if a sufficiently pre- 
cise colour curve of image 4 is obtained in the future. 


LIGHT CURVE + LENS MODEL AGE CONSTRAINTS 


In order to break the age degeneracies in the colour- 
based constraints using only data available today, we 
need to use some information about the relative bright- 
nesses of the AT 2016jka images. For this step we can 
no longer be independent of the lens models, as we 
must use the lens-model-predicted magnification values 
to de-magnify the observed photometry for comparison 
to SN models (note, however, that we again do not use 
any time delay information from the lens model for this 
method). 

For the five lens models (A-E) described above, we 
correct the observed flux density of each image (in both 
F105W and F160W bands) using the predicted lens- 
ing magnification (u). Next we employ SNTD’s “series” 
method, as it is most effective for sparse sampling, to 
attempt a reconstruction of the intrinsic SN light curve 
[59]. Once again the age of each image is constrained 
simultaneously, while also varying the SN model param- 
eters. At this stage we adopt weak priors on the intrinsic 
Type Ia SN luminosity [68] and Type Ia SN colour [69] 
to help break degeneracies in the light curve model. Ex- 
tended Data Fig. 6 shows the resulting light-curve-based 
constraints on the age of each SN image. 

As our final method to incorporate the colour and 
brightness information together, we use the colour- 
based posterior probability distributions (Extended 
Data Fig. 4) as priors for the light-curve based con- 
straints. In this approach, we must use only a single 
photometric band for the light curve constraint, so that 
we are not “double-counting” the colour information by 
simultaneously fitting to two bands together. We adopt 
F160W as the single band, since it is close to the rest- 
frame V band, where the SALT2 Type Ia SN model is 
very well constrained. The joint posterior distributions 
from this method are shown in Extended Data Fig. 7. 
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Data Availability Statement: All HST images used in this work are available from the Mikulski Archive for Space 
Telescopes (mast.stsci.edu). HST data from 2016 when the transient was active are from the program HST-GO-14496 
(archive.stsci.edu/proposal_ search.php?id=14496&mission=hst). Data collected in 2019 are from the REQUIEM 
program, HST-GO-15663 (archive.stsci.edu/proposal_search.php?id=15663&mission=hst). All VLT MUSE spectro- 
scopic data used in this work are available from the ESO Archive Science Portal (archive.eso.org). The MUSE datasets 
can be found at archive.eso.org/dataset /ADP.2019-10-07T18:14:24.762, archive.eso.org/dataset / ADP.2019-10-07T 18: 
14:24.751, and archive.eso.org/dataset /ADP.2019-10-07T18:14:24.776. All derived data supporting the findings of this 
study (photometry, lens model inputs, etc.) are available within the paper and its supplementary information files. 


Code Availability Statement: All software tools used in the analysis are publicly available, as indicated in the 
text. The software used for figure creation, including input data files, can be downloaded from github.com/gbrammer / 
mrg0138 supernova. 


Correspondence and requests for materials may be addressed to S.R. (srodney@sc.edu) and G.B. 
(gabriel.brammer@nbi.ku.dk). 
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Supplementary Material 


Telescope Instrument Band UT Date MJD Exp. Time [s] 
Spitzer IRAC 3.6 pm 2016-03-15 03:44:04 57462.156 212 
Spitzer IRAC 4.5 pm 2016-03-15 03:44:04 57462.156 241 

HST ACS/WFC F555W 2016-06-03 21:50:43 57542.910 5214 
HST WFC3/IR F160W 2016-07-18 23:14:50 57587.969* 1611 
HST WFC3/IR F105W 2016-07-19 00:43:47 57588.030* 3611 
Spitzer IRAC 3.6 pm 2016-10-13 14:35:13 57674.608 468 
Spitzer IRAC 4.5 pm 2016-10-13 14:35:13 57674.608 581 
HST WFC3/IR F110W 2019-07-13 20:53:16 58677.870 706 
AST WFC3/IR, F140W 2019-07-14 22:16:01 58678.928 353 
HST WFC3/IR F125W 2019-07-19 21:27:30 58683.894 706 
HST WFC3/UVIS F814W 2019-07-21 18:50:53 58685.785 912 
HST WFC3/UVIS F390W 2019-07-21 19:01:06 58685.792 1272 
HST WFC3/IR F140W 2019-07-21 22:42:22 58685.946 353 
VLT MUSE 0.4-0.9 um 2019-09-06 03:56:25 58732.164 2649 


Supplementary Table 1. Record of MACSJ0138 observations used in this work. Observations in which three images 
of the SN were detected are marked with an ’*’ in column five. 


Image Obs. Date (MJD) Filter Flux density [uJy] 
SN1 57588.03 F105W (Y) 0.18 + 0.02 
SN2 57588.03 F105W (Y) 2.35 + 0.02 
SN3 57588.03 F105W (Y) 0.09 + 0.02 
SN1 57587.97 F160W (H) 1.13 + 0.04 
SN2 57587.97 F160W (H) 3.57 + 0.05 
SN3 57587.97 F160W (H) 0.61 + 0.04 
SN3 58677.87 F110W (Y) 0.01 + 0.02 
SN3 58678.93 F140W (JH) 0.02 + 0.02 
SN3 58683.89 F125W (J) 0.02 + 0.03 


Supplementary Table 2. Photometry of AT 2016jka. The final three rows indicate “empty” aperture flux densities 
measured at the position of image SN3 in the 2019 HST visits. 
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ID R.A. (deg) Dec. (deg) z 

H1 24.5099018 —21.9260130 1.95 
H2 24.5132090 —21.9299032 1.95 
H3 24.5164138 —21.9303172 1.95 
H4 24.5176117 —21.9233433 1.95 
SN1 24.5151253 + —21.9306659 1.95 
SN2 24.5123198 —21.9297875 1.95 
SN3 24.5100753 —21.9273418 1.95 
3.1 24.5169659 —21.9234814 0.7663 
3.2 24.5151039 —21.9231406 0.7663 
3.3 24.5184361 —21.9264250 0.7663 
3.4 24.5146833 —21.9261020 0.7663 


Supplementary Table 3. Multiple images used as constraints in our parametric lens model. From left to right: 


image ID, right ascension, declination, spectroscopic redshift. Coordinates are in the J2000 reference frame. 
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Supplementary Table 4. Best fit model parameters for the MACS J0138 mass distribution. From left to right: mass 
component, position relative to cluster center (AR.A. and ADec.), APIE shape (ellipticity and orientation), velocity dispersion, 
core and cut radius. The final row is the generic galaxy mass at the characteristic luminosity L*, which is scaled to match each 
of cluster member galaxies. Parameters in square brackets are fixed a priori in the final model (version E). 


Method Data Lens info Priors | p(Ia) p(l)  p(Ib/c) 

a. Host color-mag Host galaxy rest-frame Hhost - 0.75 0.19 0.06 
Mx, B- K 

b. Host stellar pop. Host galaxy mass & star Hhost - 0.62 0.27 0.09 
formation rate 

c. SN color-mag SN F105W-F160W color, HSN b 0.95 0.01 0.04 
MF160W 

d. SN light curve F105W and F160WŴ SN usn, Atsn b 0.94 0.06 <0.01 


light curves 


Supplementary Table 5. SN classification probabilities for AT 2016jka. “Lens info” indicates the lensing information 
used to interpret or derive the observational data: fnost and psn are the magnifications of the host galaxy MRGO0138 and the 
SN, respectively; Atsn refers to the time delays between SN images 1, 2 and 3. In all cases the preferred LENSTOOL model E is 
used. “Priors” indicates the host galaxy classification probabilities that were adopted as priors for the subsequent classification 


using SN data. 
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SUPPLEMENTARY NOTE: COMPARISON TO PREVIOUS LENS MODELING 


Ref [7] (hereafter N18) provides a detailed analysis of this cluster, including sophisticated modeling of the lens. 
The primary modeling in N18 uses a custom lens model code [70] that traces the source Sersic profile(s) through to 
the image plane and fits the HST images directly. In addition, N18 also made a set of LENSTOOL models to estimate 
the uncertainties. Those models included both source and image plane optimization, and considered both generalized 
Navarro-Frenk-White (gNFW) [71] and dual pseudo isothermal elliptical (APIE) density profiles for the cluster. For 
this work, we used only LENSTOOL, only dPIE profiles, and considered only image-plane optimisation because it deals 
best with actual positional uncertainties. We did not incorporate pixel-level flux data, because the LENSTOOL code can 
only use observed fluxes when doing source plane optimisation (but see further discussion of this in the Supplementary 
Note: Host Image Morphology Comparison below). 

Supplementary Table 6 shows the magnification values predicted by our lens model E in comparison to the values 
reported in N18. Columns 2 and 3 give magnifications for the location of each SN image, while columns 4 and 5 
are for the host galaxy images. Column 2 repeats the SN magnifications given in the main text and Supplementary 
Table 4, which are the expectation values from the u distributions generated by the LENSTOOL MCMC sampling over 
model parameter space. Column 3 gives the “optimized” SN magnification, which is the u value returned by the single 
model instance that has the minimum y? value (note that this may be significantly different from the peak of the u 
distribution, as for image 2 in particular). Column 4 reports the minimum-y? magnification at the peak of the surface 
brightness profile for each galaxy image. Column 5 reports the ratio of the total flux to source flux of the host galaxy, 
which is effectively a flux-weighted harmonic mean of the magnification factors across each galaxy image. This last 
value is the most appropriate for comparison to the values of N18 (given in Column 6), which are computed in a similar 
way. The uncertainties from N18 reflect an estimate of systematic uncertainties, derived from lens modeling variants. 

From the last two columns of Supplementary Table 6 we see that our preferred lens model E magnifications are 
within lo of the N18 model, though our predictions are systematically lower by about 20%. This agrees with the 
assessment of systematic uncertainties from our lens model variants discussed above, which are also shown as the 
asymmetric error bars on the SN magnitudes in Fig. 2 in the main text, and in Extended Data Fig. 2. As seen in 
those figures, a larger magnification value would make the SN more consistent with the expected luminosity of a Type 
Ia SN at this redshift—though it would also make it more consistent with some CC SN light curves, likely making the 
classification somewhat more ambiguous. 

It is important to note that the all of the model variants explored here adopted dPIE distributions as the density 
profile for all lensing components. It is well-documented that the choice of the density profile can strongly affect 
the inferred magnifications and time delays in a strong lensing system. This is therefore another potential source of 
systematic uncertainty that should be explored in future lens modeling. The models of N18 explored models using 
gNFW profiles, but of course did not make explicit predictions for the magnifications or time delays at the SN locations. 
If alternate density profiles result in significant changes to the model-predicted magnfications for the SN images, that 
could in principle change the conclusions about SN classification and age constraints described here. 

We expect that new lens modeling of this cluster with alternate software and different choices of constraints will also 
be informative, and may improve on our model predictions. Both the N18 modeling and the construction of lens model 
E and our model variants are constructed blind (without knowledge of the SN magnitudes). Future lens modeling 
could incorporate the measured magnitudes as constraints, potentially yielding a more robust prediction of the time 
delay for the fourth image. We hope that the discovery of AT 2016jka will encourage such efforts. 


SUPPLEMENTARY NOTE: HOST IMAGE MORPHOLOGY COMPARISON 


As a qualitative check to evaluate the accuracy of lens model E, we also simulated the overall shape of the SN host 
galaxy using the sum of two Sersic profiles, and then propagated this through the lens model to make predictions for 
the morphology of the host arcs. The results are shown in Supplementary Fig. 1 for host images H1-H3 (1.1-1.3). Here 


Image ( HSN ) Hopt,SN Hopt,gal.pk. Hgal.avg. HN18,gal.avg. 


1 3.90.5 6.7 4.35 10.0 12.5+5.4 
2 T.43 15.2 6.9 8.3 10.33.1 
3 51 4.3 3.64 4.2 4.9+1.6 


Supplementary Table 6. Comparison of magnification predictions with the N18 lens model. 
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Supplementary Figure 1. Comparison of the SN host galaxy against a simulation using lens model E. North is 
up and East is left. Image 1.1 is shown in the top row, 1.2 in the middle, and 1.3 in the bottom. The left column (panels 
a,d,g) shows each observed image in the F160W bandpass, with a scale bar indicating 1 arcsecond. The black cross highlights 
the location of the SN image in each case. The central column (b,e,h) shows a simulated galaxy profile generated using lens 
model E, with a sum of two Sersic profiles to represent the background galaxy. The last column (c,f,i) shows the residuals 
(observed-model). Overlaid contours in the first column show the profile of the observed data (green) and the model (orange). 
The contours are logarithmically spaced and correspond to levels of 1, 2, 4, and 80 in detection, where ø is the pixel to pixel 
background level. Some residual flux is apparent in the difference images, but the overall shape of each image is well matched 
by the simulated profile and the distortions introduced by lens model E. 


Supplementary Figure 2. Comparison of the SN host galaxy image 4 against a simulation using lens model E. 
Same as Figure 1, but showing galaxy image H4 (1.4). 
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we see the model produces a very good match to the observed arcs. This is especially encouraging because the surface 
brightness distributions of the arcs were not used as inputs for the lens modeling. 

The predicted arc morphology for host galaxy image H4 (1.4) is shown in Supplementary Fig. 2. The residuals are 
more significant, showing clear asymmetric structure, which may indicate a mismatch of rotation and/or shear angle 
at the location of this image. Nevertheless, the global morphology of the arc is similar in location, size and elongation, 
which again could be taken as an encouraging indicator of the validity of magnification and time delay estimates for 
image SN4. For image 5 (not shown), lens model E predicts a very compact source at the location of the BCG, which 
is at odds with the clearly elongated radial arc of image 5 that can be seen in the HST imaging. This may well be 
due to the fact that the lens model overestimates the velocity dispersion of the galaxy (see above, and Supplementary 
Note: Future Lens Modeling). Reconstruction of image H5 is not as important, because the highly demagnified final 
transient image SN5 is not expected to be observable anyway. Nonetheless, a more accurate reconstruction of host 
galaxy image H4 (and to a lesser extent image H5), should be an important metric for future lens modeling. 


SUPPLEMENTARY NOTE: FUTURE WORK 
FUTURE LENS MODELING 


The uncertainty in the predicted time delay from the lens model presented here is +540 days (in the observer frame). 
Are there improvements to the lens modeling that could tighten this prediction? 

Our lens modeling did not make use of the velocity dispersion of the BCG, measured as 390+10 km s~! from our 
MUSE data. We did not have this measurement available prior to unblinding, and thus we have restricted ourselves 
to only consider fully blind lens models in this paper. However, we note that this measured value is very discrepant 
with the best-fit value of 700 km s~! from our preferred lens model, variant E (see Supplementary Table 4). Future 
analyses could incorporate the BCG velocity dispersion and here we speculate about the impact this would have. 

With a preliminary extension of model E we find that the time delay predictions would likely shift by <15 days, and 
the predicted magnifications for SN images 1-3 would likely be larger by as much as ~2.5x. This is within the range 
of our systematic uncertainty estimates based on lens models A-E. We therefore expect that inclusion of the BCG 
velocity dispersion will not significantly impact the SN classification or time delay conclusions, but may improve the 
accuracy and precision of the time delay estimate for image SN4. 

One may anticipate that future observations of the lensing cluster (e.g. with JWST) will provide more lens modeling 
constraints, potentially including the discovery of new multiply-imaged systems and new redshift measurements. These 
would substantially improve the lens model time delay predictions. Application of different lens modeling approaches 
could also provide some added confidence that there are not systematic biases in the lens model time delay prediction. 


FUTURE OBSERVATIONS 


Though the +540 day time delay uncertainty is small relative to the baseline of > 7000 days, it nevertheless represents 
a long period over which the field would need to be monitored for the reappearance. Is it feasible to expect future 
observations to catch the reappearance of AT 2016jka close to peak brightness? Even if the time delay uncertainty 
remains on the order of +1 year, it would still be reasonable to execute a follow up campaign with a cadence of 
approximately 2 months. Since the SN is at a redshift of z = 1.95, a span of, say, 50 days in the observer frame 
is 17 days in the rest frame, comparable to the rise-time of a Type Ia SN. Thus, it is reasonable to expect that a 
relatively inexpensive monitoring campaign would be able to catch the return appearance of AT 2016jka at or before 
peak brightness, ensuring a well-sampled light curve for the final SN image. 

If future follow-up observations are successful in capturing the full light curve of the final fourth image, then future 
lens modeling could also incorporate measured SN magnifications as constraints. Measured magnifications from lensed 
Type Ia SNe have been used to test lens models at both the cluster and galaxy scale [72; 73; 74]. The addition of 
astrometric constraints for the fourth image could also significantly improve the time delay predictions for a lens model 
[75]|—and this could be done even with a fully blinded analysis. Measurement of the magnification for a lensed Type Ia 
SN can be done without adopting strong priors from a cosmological model [76], meaning that one can avoid a circular 
constraint when the AT 2016jka time delay is then used for cosmology. 


FUTURE DISCOVERIES 


In addition to follow-up observations of AT 2016jka, we may also hope for more discoveries of similar cluster-lensed 
SNe with long time delays. A primary motivation for pursuing such events is that they can be a relatively low-cost 
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tool for time delay cosmography. As AT 2016jka shows, when the time delay is longer than a few years, the time delay 
measurement can be anchored at either end by just a few epochs of imaging. If similar events are detected while the 
SN is still observable, one could collect a well-sampled light curve for an early and bright image using ground-based 
telescopes. After waiting through the decade-long delay, the SN’s reappearance can be captured with a relatively low- 
cost monitoring campaign. A full light curve of the final image would not be needed. For example, with AT 2016jka 
even if the time delay is measured to only +150 days, that would be a 2% time-delay measurement, meaning it is 
lens-model-limited for cosmological inferences. 

The expected rate for such events is still highly uncertain, and published rate estimates to date can only be taken 
as extreme lower limits for the expected yield from future sky surveys |77; 46; 78; 79]. All of these past analyses have 
been limited to < 5 well-studied galaxy clusters. Furthermore, they have only examined the set of already known 
multiply-imaged galaxies, and have explicitly predicted only the rate of events that would have a time delay of <5 
years. With these caveats, the predicted lower limits are of order 1 SN detection per year per cluster, for a deep survey 
with a detection limit of 27 AB mag [46]. At the 50 limits of the Rubin Observatory (i ~ 23.4 AB mag), the lower 
limit on that rate is reduced by about a factor of ten [79]. 

It is treacherous to extend these estimates to the larger population of all galaxy clusters that will be regularly observed 
by future wide-field surveys. Nevertheless, let us make a crude extrapolation to motivate future work. Consider the 
1-year 2000 deg? High Latitude Survey (HLS) from the Roman Space Telescope [19; 80], and let us conservatively 
apply the rate of ~1 SN yr! cluster! to only the ~10 most massive clusters in the HLS area. This still predicts 
at least 10 cluster-lensed SN detections, which is comparable to the few dozen galaxy-lensed SNe expected from the 
Roman SN cosmology survey [81]. Similarly, if we apply the 10x lower rate for the Rubin Observatory to the most 
massive clusters in the LSST survey area, we would anticipate at least ~10 detections over the 10-year survey. This 
discovery rate from wide-field surveys could be enhanced with dedicated ground-based cluster surveys. [82; 83; 77] 

We hope that the discovery of AT 2016jka will motivate an improvement over this very rough estimation of future 
rates. This would require a more complete census of lensing clusters, along with lens models to predict magnifications 
and time delays, and measurements of star formation and stellar mass in lensed galaxies to predict the SN explosion 
rates. 


REFERENCES 


[78]Petrushevska, T. et al. Searching for supernovae in the 
Dark Matter Profiles from Galaxies to Clusters: Bridging 
the Gap with Group-scale Lenses. ApJ 814, 26 (2015). 


multiply-imaged galaxies behind the gravitational 
telescope A370. A&A 614, A103 (2018). 


[79]Petrushevska, T. Strongly lensed supernovae in 


278, 488 (1996) well-studied galaxy clusters with the Vera C. Rubin 


Observatory. Symmetry 12, 1966 (2020). 


cluster mass models: MNRAS 440,242 (2014), [80]Troxel, M. A. et al. A synthetic Roman Space Telescope 


. . High-Latitude Imaging Survey: simulation suite and the 
supernova magnified by the frontier fields galaxy cluster 


abell 2744. ApJ 811, 70 (2015). impact of wavefront errors on weak gravitational lensing. 


MNRAS 501, 2044-2070 (2021). 


constraints from the first resolved strongly lensed Type Ia [81]Pierel, J. D. R. et al. Projected Cosmological constraints 


supernova iPTFl6geu. MNRAS 491, 2639-2654 (2020). 

75|Birrer, S. & Treu, T. Astrometric requirements for strong 
lensing time-delay cosmography. MNRAS 489, 
2097-2103 (2019). 

76|Patel, B. et al. Three gravitationally lensed supernovae 
behind CLASH galaxy clusters. ApJ 786, 9 (2014). 

77|Riehm, T. et al. Near-IR search for lensed supernovae 
behind galaxy clusters. III. Implications for cluster 
modeling and cosmology. A&A 536, A94 (2011). 


from strongly lensed supernovae with the Roman Space 
Telescope. ApJ 908, 190 (2021). 

[82]Stanishev, V. et al. Near-IR search for lensed supernovae 
behind galaxy clusters. I. Observations and transient 
detection efficiency. A&A 507, 61-69 (2009). 

[83]Goobar, A. et al. Near-IR search for lensed supernovae 
behind galaxy clusters. II. First detection and future 
prospects. A&A 507, 71-83 (2009). 


